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ABSTRACT 

Aims. We study the changes in the properties of turbulence driven by the magnetorotational instability in a shearing box, as the 
computational domain size in the radial direction is varied relative to the height 

Methods. We perform 3D simulations in the shearing box approximation, with a net magnetic flux, and we consider computational 
domains with different aspect ratios 

Results. We find that in boxes of aspect ratio unity the transport of angular momentum is strongly intermittent and dominated by 
channel solutions in agreement with previous work. In contrast, in boxes with larger aspect ratio, the channel solutions and the 
associated intermittent behavior disappear. 

Conclusions. There is strong evidence that, as the aspect ratio becomes larger, the characteristics of the solution become aspect ratio 
independent. We conclude that shearing box calculations with aspect ratio unity or near unity may introduce spurious effects. 
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1. Introduction 

Understanding the process of angular momentum transport is 
one of the fundamental issues in the physics of accretion disks. 
Turbulence was reco gnized as the main s ource of the required 
enhanced transport dShakura & Sunvaevl Il973|). but its driv- 
ing mechanism remained unclear until Balbus & Hawlevl d!991l) 
proposed the magnerotational instability (MRI). Much of what 
is presently known about MRI driven turbulence has been ob- 
tained in the shearing box approximation in whic h only a small 
perio d ic patch of the di s k is s i mulated(see e. g . lHawlev et all 
1 19951 ISano & Inutsukal l200lt ISano & Stond 120021: iBalbusl 
bpoalSano et al.Ll2004tlTurner et al.ll2003l:lLesur & Longarettil 
2007). The advantage of this local approach, as opposed to a 
global disk simulation, is the possibility of reaching much higher 
resolutions at the same computational cost. The shearing box 
approach is of course meaningful only if it captures the char- 
acteristics of MRI driven turbulence, and of the related angular 
momentum transport in a full disk. A critical discussion of the 
validity of th e shearing box approxima tion can be found in a re- 
cent work by Regev & Umurhan (2007) 

Ideally one should compare local and glob al results ; how- 
ever available global disk sim ulations (see e.g. iHawlevi 120011: 
lHawl ev. Balbus & Stone], 1200 lb have relatively low resolution 
while simulations with adequately high resolution are still be- 
yond present capabilities. Given the present resources, one 
should minimally check the self consistency of the shearing box 
results. The shearing box equations are obtained as a for mal lo- 
cal ex pansion in the limit of large radii and small gap dHill L 
Il878l) . but there is no guarantee that the solutions thus obtained 
satisfy the same locality conditions. This should be verified a 
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posteriori, by checking, for example, that the properties of the 
solutions do not depend on the size of the computational domain. 

It useful to distinguish two related issues. One concerns the 
size of the computational domain relative to some characteristic 
length scale of the problem determined by the physical param- 
eters. For instance, in the MRI case, the vertical wavelength of 
the mode of maximum growth rate depends on the strength of the 
applied uniform magnetic field. For a given field strength then, 
the vertical size of the computational domain is chosen to con- 
tain some multiple of that wavelength. In this case, varying the 
field strength is in some sense equivalent to varying the (vertical) 
box size. A related issue concerns the aspect ratio, i.e the size of 
the computational domain in the radial and azimuthal directions 
relative to the vertical size. The modes of maximum growth rate 
do not offer much guidance in this respect since they only vary in 
the vertical. Clearly a very thin box may introduce spurious ef- 
fects, a very wide box rapidly becomes computationally expen- 
sive. A typical compromise is to adopt boxes with aspect ratio 
between the radial and vertical direction of unity. Presumably, 
the basis for this choice is that once nonlinear processes satu- 
rate the exponential growth, the typical resulting structures will 
have a characteristic size in the radial direction comparable, or 
in any case, not much larger than that in the vertical. This being 
the case, a computational domain of unit aspect ratio is adequate 
to capture the properties of the solution. Strangely, this point has 
not received careful consideration and there has been no system- 
atic study of the dependence of the solutions on the aspect ra- 
tio. Furthermore many studies in computational domains of unit 
aspect ratio have shown intermittent behavior associated with 
the formation and subsequent disruption of "chan nel" solutions 
dHawlev & Balbuslll992l:ISano & Inutsukall2001h . These are the 
nonlinear analogues of the exponentially growing linear modes. 
Recently a number of authors have addressed the effect of the 
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presence or absence of channel solutions in different contexts, 
for example, by introd ucing vertical gravity and stratification 
( Cop pi & Kevesl 120031) or changing the radial boundary con- 
ditions (|Liu. G oodman &Ji, 2006; Umurhan, Regev & M enoul 
l2007tlUmurhan. Menou & Regevil2007l) . Here we note that, in- 
terestingly, in the nonlinear regime, the observed channel solu- 
tion is not necessarily the one that corresponds to the linear mode 
of maximum growth rate. In fact, often, the observed channel 
solution is the one whose vertical extent fills the computational 
domain once. In this case, it is not a priori obvious what as- 
pect ratio should be used. In this paper we address these issues 
by carrying out a systematic study of shearing box results as a 
function of aspect ratio. In particular, we want to see if the re- 
sults observed in boxes of unit aspect ratio are representative of 
more extended systems, or if they display peculiarities induced 
by an overly constrained geometry. 
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Fig. 2. Probability distribution functions of the Maxwell stresses for 
the same cases as in Figure fJJ The three curves correspond to the three 
values of the aspect ratio: L = 1 — solid curve, L = 4 — dashed curve, 
and L = 8 — dot-dashed. 



2. Formulation 

We perform a series of 3D, compressible, isothermal numerical 
simulations in the shearing box approximati on (for a detailed d e- 
scription of the shearing box model, see lHawlev et al.1 SmEti ). 
The Cartesian coordinates x, y, z are defined around a fiducial 
point comoving with angular velocity Q. and correspond, respec- 
tively, to the radial, azimuthal and vertical directions. An equi- 
librium solution is given by a uniform shear flow, v = -qQ.xe y , 
with q — 3/2 for a Keplerian disk, threaded by a uniform ver- 
tical magnetic field, B = Boe : , in which the density p and the 
pressure p are constant (p - 1, p = c 2 s p). We have assumed, as 
it is done by most of the literature on the subject, the simplest 
case with no vertical stratific ation and gravity. This ma y limit 
the applicability to real disks dRegev & Umurhanl 120071) . but it 
is adequate for our purpose of studying the aspect ratio depen- 
dence. This equilibrium is initially perturbed by small amplitude 
spatially uncorrelated azimuthal velocity perturbations. The box 
has size L x 4 x 1 in the x, y and z directions, respectively, with 
L varying from 1 to 8. 

We take 1 /Q. as the unit of time (note however that in the 
plots time is in unit of the rotation time 2w/Q), and in all the 
simulations the sound speed c s and the plasma B = p gas /Pmag 
are fixed with values of 4.56 and 10 4 , respectively. With these 
choices the fastest growing MRI mode has a vertical wavelength 
close to 1/3. 

Computations are carried at both low and high resolutions 
(32 and 128 zones per unit length, respectively) on equally 
spaced grids. The MHD equations are solved in conservative 
form using the isothermal MHD modu l e ava ilable in the PLUTO 
code dMignone et all l2007t iMignond 120071) . The latter is a fi- 
nite volume, Riemann solver based code in which the evolu- 
tion of the magneti c field is carried out us ing the constrained 
transport method of Balsara & Spicer ( 1999). In the present for- 
mulation we do not include explicit viscosity or magnetic dif- 
fusivity, and rely solely on numerical dissipation to limit the 
potential unbounded growth of the solutions (for a theoretical 
discussion of numerical dis sipati o n in Riemann solver based 
schemes see, for example, iTorol dl999t) : for a discussion of 
the role of numerical dissipation i n turbulent simulations see 
iGrinstein, Margolin. & Riderld2007l) ). 

3. Results 

Our main objective is to study how the properties of the turbu- 
lent solutions vary as the (radial) aspect ratio L is varied from 1, 



to 4, to 8. We recall that in simulations with aspect ratio unity, 
the angular momentum transport shows an intermittent behav- 
ior with episodes of high transport. These correspond to states, 
loosely referred to as "channel" solutions, in which the veloc- 
ity and magnetic perturbations are highly spatially correlated. 
Strictly speaking, the channel solutions are exponentially grow- 
ing exact solutions of the incompressible shearing box equa- 
tions in which the velocity and magnetic field perturbations de- 
pend sinusoidally on z, have directions at right angle to each 
other, and a growth rate determined by the ang le they make 
to the azimuthal direction dGoodman & Xul 1 19941) . More com- 
monly the term "channel solutions" is used to describe a behav- 
ior with similar prope rties to that observed i n the exact solu- 
tions. It was shown by Goodman & Xul d 19941) that the channel 
solutions, as their amplitudes grow, become unstable to para- 
sitic instabilities. It is now believed that the formation of near 
channel solutions and their subsequent disruption by the para- 
sitic in stabilities are the basis f or the observed intermittent be- 
havior dSano & Inutsukall2001l) . It should be noted that the ver- 
tical wavenumber determines the angle y c between the direction 
of the magnetic perturbations and the azimuthal direction, and 
hence the growth rate. Although, in principle, channel solutions 
exist with any wavenumber, the most commonly observed ones 
have a vertical extent equal to the vertical box size, for which the 
angle y c w 23°. We shall use these results presently. 

We start our discussion with the low resolution results and 
compare three cases with different aspect ratio. Fig.Q]shows the 
the volume averaged Maxwell stresses (w xy ) = -(B x B y ) as a 
function of time. In what follows we concentrate on the Maxwell 
stresses because they give the largest contribution to the total 
stresses exceeding the Reynolds stresses on average by a fac- 
tor of about 5. We have normalized the stresses to the pressure, 
thu s the values shown in the fi gure correspond to the a parame- 
ter dShakura & Sunvaevlll973l) . The left panel, corresponding to 
L = 1, shows the characteristic intermittent behavior described 
above. A striking difference appears if we compare this case with 
the other two having larger aspect ratios in which the spikes are 
absent. A more quantitative view of the differences between the 
three cases is afforded by their respective probability distribution 
functions shown in Fig [2] The case L = 1 (solid curve) is quite 
distinct from the other two with a long tail corresponding to the 
episodes of enhanced transport. Careful examination shows that 
a small tail is still present in the case with L = 4 (dashed curve), 
and absent in the case with L = 8 (dot-dashed curve). However, 
we see that the differences between the last two cases are much 
smaller than the differences with the L = 1 case, indicating that 
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Fig. 1. Time histories of the Maxwell stresses averaged over the computational box and normalized to the pressure for the low resolution simula- 
tions. The three panels correspond to cases with L = 1,4, and 8. Note that time is in unit of the rotation time. 



we have almost reached a converged behavior. Because of the 
presence of the peaks, the time-averaged value of a for L = 1 
exceeds the other two cases by approximately a factor of 2. It is 
natural to assume that the distinctive behavior of the case with 
L = 1 is due to the presence of channel solutions that are other- 
wise absent in the simulations with larger aspect ratios. This can 
be verified by looking for typical imprints of the channel solu- 
tions such as, for instance, a high correlation coefficient between 
directional components of the magnetic filed, or a specific rela- 
tionship between the vertical wavenumber and the angle between 
the magnetic field and the azimuthal direction. In order to define 
these quantities it is useful to introduce the idea of a scatter dia- 
gram whereby for each gridpoint the value of B y is plotted as a 
function of B x , Fig [4] For a pure channel solution this procedure 
would yield a straight line through the origin with slope tan y c . In 
the general case, the points will not lie on a straight line; still, a 
least-square straight line fit through the points can be used to de- 
fine the average angle, and the spread about this line would give 
a measure of the correlation coefficient R. Equivalently, these 
quantities can be defined by 



tan y = — 



( BxBy) 



and R 



(B X ByY 



(1) 



The time histories of these quantities for different aspect ratios 
are shown in Fig. [3] clearly indicating the presence of channel 
solutions in the case with aspect ratio unity and their absence 
in the other cases. Comparison between Figures Q] and [3] also 
confirms that the spikes with high transport correspond, indeed, 
to the formation of channel solutions. 

Differences between the various cases can be further illus- 
trated by scatter plots of the distribution of B y as a function of B x 
like those presented in Fig. [4] The two panels on the left show 
two examples of distributions for the case L = 1, one corre- 
sponding to a maximum of (w xv ) and the other corresponding to 
a minimum. This confirms that in going from the minimum to the 
maximum the magnetic field fluctuations increase their intensity 
and the x and y components tend to become more correlated. 
The velocity fluctuations show a similar behavior, with almost 
no correlation in the minimum state. In contrast, for L = 4, there 
are no significant variations in the distributions during the evo- 
lution, see Fig. |4] where we display only one representative ex- 
ample that shows a similarity with the minimum state of the case 
with L — 1 . The figure also makes apparent why the Reynolds 
stresses are significantly lower than the Maxwell stresses. The 
difference partly arises because the intensity of the fluctuations 
is smaller, but to a somewhat larger extent because the veloc- 
ity fluctuations, except during a channel solution phase, remain 
uncorrelated. 

Another property of the distributions not entirely evident 
from the figure is that, while for L = 4 and at the minimum of the 
L = 1 case the distributions peak at zero value, at the maximum 





Fig. 4. Scatter diagrams of horizontal magnetic fluctuations (dark 
tones) and velocities (lighter tones). The first panel corresponds to the 
case with L = 1 at an instant near a maximum in the stress. The second 
panel corresponds to the same case near a minimum. The third panel 
corresponds to a representative instant for the case with L = 4. 



of the L — 1 case there are two peaks around the largest values 
of the fluctuations. This is clear in Fig. [5] where we show the 
probability distribution function of the azimuthal magnetic field 
component B y . The panels show the behavior at an instant near a 
maximum in the stress (solid line) and a minimum (dashed line) 
for two different aspect ratios. In the left panel, we also show 
(dot-dashed line) the corresponding distribution for the channel 
solution. The presence of the double peak distribution in the case 
with L — 1 (left panel) can be considered again as an indica- 
tion of the presence of the channel solution around the maxima 
in the stress, as it can be seen by com paring the solid and the 
dot-dashed lines. Recently, the work of iFromang & Papaloizoul 
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Fig. 5. Probability distribution functions for the azimuthal component 
of magnetic field B y . The left panel refers to the case with L = 1. Solid 
and dashed lines correspond to the first two panels in Fig. [4] while the 
dot-dashed line corresponds to the exact channel solution. The right 
panel corresponds to the third panel in Fig. [4] 



(2007) has raised some concerns about the dependence of nu- 
merical studies of MRI turbulence on resolution. Accordingly, 
we have repeated the simulations with L — 1 , 4 at a higher reso- 
lution of 128 gridpoints per unit length. We find the results to be 
qualitatively the same. For example, Fig.|6]shows the time histo- 
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Fig. 3. Time histories of the average slope tan y (upper row) and correlation coefficient R (lower row) for the same cases as in Figure Q] The 
straight horizontal lines indicate the theoretical value of tan y c for the exact channel solution with wavelength unity. Note that time is in unit of the 
rotation time. 
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Fig. 6. Time histories of the Maxwell stresses averaged over the com- 
putational box and normalized to the pressure for the high resolution 
simulations. The two panels correspond to cases with L = 1, and 4. 
Note that time is in unit of the rotation time. 



ries of (w xy ) for the higher resolution simulations. The behavior 
is strikingly similar to that of Fig. [1] with a strong presence of 
the channel solutions in the L = 1 case and a complete absence 
of such solutions in the runs with larger aspect ratio. We have 
also checked other indicators, like the probability distribution 
functions and the scatter diagrams and found them to be qual- 
itatively the same as in the cases described above. We should 
note however that quantitatively we do observe a dependence 
on resolution; in particular the average value of the Maxwell 
stress es increases with increasing resolution (cf. the zero flux 
case Fromang & Papaloizou (2007)). This follows mainly from 
an increase in the peak values of magnetic fluctuations, which 
is not surprising since the effective magnetic Reynolds number 
associated with numerical dissipation increases with increasing 
resolution. 



4. Conclusion 

Channel solutions represent an highly correlated state for which 
the transport is very efficient and are a characteristic feature 
that strongly affects the behavior of the angular momentum 
transport in most of the shearing box simulations so far pre- 



sented in the literature. There is however the question whether 
the dominance of this state could not be due to the constraints 
imposed by the typical shearing box approach. Recently, sev- 
eral authors have addressed this questio n considering th e effect s 
of vertical gravity and stratification (Coppi & Keves, 2003), 
and of different boundary conditio n s dLiu. G oodman &Ji, 2006; 
Umurhan, Menou & Regev, 2007; Umurhan, Regev & Menou, 
2007). The general result is that the channel solution is not ro- 
bust and disappears whenever one introduces some modification 
to the typical shearing box approach and, as a consequence, the 
angular momentum transport is less efficient. 

In this paper we have focused our analysis on the effects of 
a change of the box aspect ratio and the results above show two 
different behaviors depending on it. For aspect ratios close to 
unity the solutions are strongly affected by the channel solutions, 
with frequent episodes of high correlation and efficient trans- 
port. For larger aspect ratios the channel solutions disappear, the 
system remains in the state with lower correlation and the aver- 
age transport of angular momentum is correspondingly reduced. 
Further increase in aspect ratio does not lead to any significant 
changes. Thus we conclude that the less correlated state is more 
likely to be representative of the extended system. 

It is interesting to speculate why the channel solutions disap- 
pear from boxes with larger aspect ratio. One possibility is to as- 
sume that nontrivial correlations are necessary to form the chan- 
nel solutions and that these cannot develop faster than the sound 
crossing time. Accordingly, in a large box it takes longer to form 
the channel solutions; if the rate at which they are destroyed is 
fixed in a sufficiently large box they may never form. We do not 
believe this is the correct explanation, since we have found sig- 
nificantly the same behavior in simulations with a sound speed 
ten times larger. Another, more likely explanation, is that the sec- 
ondary instabilities responsible for the destruction of the channel 
solutions are to some extent suppressed at small aspect ratios, 
but can achieve their maximal growth rate once the aspect ratio 
is large enough. In supp ort of this idea, we n ote that the parasitic 
instabilities studied bv lGoodman & Xul (1 19941) in general have a 
horizontal size larger than the vertical size of the channel solu- 
tion on which they grow. In a box of unit aspect ratio containing 
a single channel they may simply be stable. 

Whatever the reason for the differences, clearly, the obvious 
conclusion is that studies of turbulence driven by the MRI with 
net flux should be conducted in shearing boxes sufficiently large 
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to allow the solution to develop naturally. The exact domain size 
depends on the particular problem and the quantities under con- 
sideration. However, it seem clear that boxes with aspect ratio 
close to unity over-emphasize the role of the channel solutions 
and may lead to misleading results. 

On the other hand our discussion does not put in doubt the 
fact that MRI may be adequate to produce the turbulence nec- 
essary to support the required "viscosity" to transport angular 
momentum. However, the uncertainties discussed by many au- 
thors about the applicability of the shearing box approximation 
for the astrophysical problem of accretion disks suggest that the 
extension of the shearing box results to the full disk has to be 
tested on global simulations with sufficient resolution, that may 
soon be feasible. 
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